{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nNetwork group process\n=====================\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 2\nimport re\nfrom datetime import datetime\n\nfrom matplotlib import pyplot as plt\nfrom matplotlib.animation import FuncAnimation\nfrom matplotlib.colors import ListedColormap\nfrom numba import njit\nfrom numpy import arange, array, empty, ones\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.generic import flatten_line_matrix\nfrom py_eddy_tracker.observations.network import Network\nfrom py_eddy_tracker.observations.observation import EddiesObservations"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class VideoAnimation(FuncAnimation):\n    def _repr_html_(self, *args, **kwargs):\n        \"\"\"To get video in html and have a player\"\"\"\n        content = self.to_html5_video()\n        return re.sub(\n            r'width=\"[0-9]*\"\\sheight=\"[0-9]*\"', 'width=\"100%\" height=\"100%\"', content\n        )\n\n    def save(self, *args, **kwargs):\n        if args[0].endswith(\"gif\"):\n            # In this case gif is used to create thumbnail which is not used but consume same time than video\n            # So we create an empty file, to save time\n            with open(args[0], \"w\") as _:\n                pass\n            return\n        return super().save(*args, **kwargs)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "NETWORK_GROUPS = list()\n\n\n@njit(cache=True)\ndef apply_replace(x, x0, x1):\n    nb = x.shape[0]\n    for i in range(nb):\n        if x[i] == x0:\n            x[i] = x1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Modified class to catch group process at each step in order to illustrate processing\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class MyNetwork(Network):\n    def get_group_array(self, results, nb_obs):\n        \"\"\"With a loop on all pair of index, we will label each obs with a group\n        number\n        \"\"\"\n        nb_obs = array(nb_obs, dtype=\"u4\")\n        day_start = nb_obs.cumsum() - nb_obs\n        gr = empty(nb_obs.sum(), dtype=\"u4\")\n        gr[:] = self.NOGROUP\n\n        id_free = 1\n        for i, j, ii, ij in results:\n            gr_i = gr[slice(day_start[i], day_start[i] + nb_obs[i])]\n            gr_j = gr[slice(day_start[j], day_start[j] + nb_obs[j])]\n            # obs with no groups\n            m = (gr_i[ii] == self.NOGROUP) * (gr_j[ij] == self.NOGROUP)\n            nb_new = m.sum()\n            gr_i[ii[m]] = gr_j[ij[m]] = arange(id_free, id_free + nb_new)\n            id_free += nb_new\n            # associate obs with no group with obs with group\n            m = (gr_i[ii] != self.NOGROUP) * (gr_j[ij] == self.NOGROUP)\n            gr_j[ij[m]] = gr_i[ii[m]]\n            m = (gr_i[ii] == self.NOGROUP) * (gr_j[ij] != self.NOGROUP)\n            gr_i[ii[m]] = gr_j[ij[m]]\n            # case where 2 obs have a different group\n            m = gr_i[ii] != gr_j[ij]\n            if m.any():\n                # Merge of group, ref over etu\n                for i_, j_ in zip(ii[m], ij[m]):\n                    g0, g1 = gr_i[i_], gr_j[j_]\n                    apply_replace(gr, g0, g1)\n            NETWORK_GROUPS.append((i, j, gr.copy()))\n        return gr"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Movie period\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t0 = (datetime(2005, 5, 1) - datetime(1950, 1, 1)).days\nt1 = (datetime(2005, 6, 1) - datetime(1950, 1, 1)).days"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Get data from period and area\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "e = EddiesObservations.load_file(data.get_demo_path(\"network_med.nc\"))\ne = e.extract_with_mask((e.time >= t0) * (e.time < t1)).extract_with_area(\n    dict(llcrnrlon=25, urcrnrlon=35, llcrnrlat=31, urcrnrlat=37.5)\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reproduce individual daily identification(for demonstration)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "EDDIES_BY_DAYS = list()\nfor i, b0, b1 in e.iter_on(\"time\"):\n    EDDIES_BY_DAYS.append(e.index(i))\n# need for display\ne = EddiesObservations.concatenate(EDDIES_BY_DAYS)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run network building group to intercept every step\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = MyNetwork.from_eddiesobservations(EDDIES_BY_DAYS, window=7)\n_ = n.group_observations(minimal_area=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def update(frame):\n    i_current, i_match, gr = NETWORK_GROUPS[frame]\n    current = EDDIES_BY_DAYS[i_current]\n    x = flatten_line_matrix(current.contour_lon_e)\n    y = flatten_line_matrix(current.contour_lat_e)\n    current_contour.set_data(x, y)\n    match = EDDIES_BY_DAYS[i_match]\n    x = flatten_line_matrix(match.contour_lon_e)\n    y = flatten_line_matrix(match.contour_lat_e)\n    matched_contour.set_data(x, y)\n    groups.set_array(gr)\n    txt.set_text(f\"Day {i_current} match with day {i_match}\")\n    s = 80 * ones(gr.shape)\n    s[gr == 0] = 4\n    groups.set_sizes(s)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Anim\n----\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(16, 9), dpi=50)\nax = fig.add_axes([0, 0, 1, 1])\nax.set_aspect(\"equal\"), ax.grid(), ax.set_xlim(26, 34), ax.set_ylim(31, 35.5)\ncmap = ListedColormap([\"gray\", *e.COLORS[:-1]], name=\"from_list\", N=30)\nkw_s = dict(cmap=cmap, vmin=0, vmax=30)\ngroups = ax.scatter(e.lon, e.lat, c=NETWORK_GROUPS[0][2], **kw_s)\ncurrent_contour = ax.plot([], [], \"k\", lw=2, label=\"Current contour\")[0]\nmatched_contour = ax.plot([], [], \"r\", lw=1, ls=\"--\", label=\"Candidate contour\")[0]\ntxt = ax.text(29, 35, \"\", fontsize=25)\nax.legend(fontsize=25)\nani = VideoAnimation(fig, update, frames=len(NETWORK_GROUPS), interval=220)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Final Result\n------------\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(16, 9))\nax = fig.add_axes([0, 0, 1, 1])\nax.set_aspect(\"equal\"), ax.grid(), ax.set_xlim(26, 34), ax.set_ylim(31, 35.5)\n_ = ax.scatter(e.lon, e.lat, c=NETWORK_GROUPS[-1][2], **kw_s)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}